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ABSTRACT 

We measure the mass function of dark matter halos in a large set of collisionless cosmological simulations 
of flat ACDM cosmology and investigate its evolution at z < 2. Halos are identified as isolated density peaks, 
and their masses are measured within a series of radii enclosing specific overdensities. We argue that these 
spherical overdensity masses are more directly linked to cluster observables than masses measured using the 
friends-of-friends algorithm (FOF), and are therefore preferable for accurate forecasts of halo abundances. Our 
simulation set allows us to calibrate the mass function at z = for virial masses in the range 10 11 h~ l M Q 
< M < 10 15 h~ l Mq to < 5%. We derive fitting functions for the halo mass function in this mass range for 
a wide range of overdensities, both at z = and earlier epochs. In addition to these formulae, which improve 
on previous approximations by 10-20%, our main finding is that the mass function cannot be represented 
by a universal fitting function at this level of accuracy. The amplitude of the "universal" function decreases 
monotonically by « 20-50%, depending on the mass definition, from z = to 2.5. We also find evidence for 
redshift evolution in the overall shape of the mass function. 

Subject headings: cosmology:theory — dark mattenhalos — methods:numerical — large scale structure of the 



1. INTRODUCTION 

Galaxy clusters are observable out to high redshift (z < 1— 
2), making them a powerful probe of cosmology. The large 
numbers and high concentration of early type galaxies make 
clusters bright in optical surveys, and the high intracluster gas 
temperatures and densities make them detectable in X-ray and 
through the Sunyaev-Zel'dovich (SZ) effect. The evolution 
of their abundance and clustering as a function of mass is 
sensitive to the power spectrum normalization, matter con- 
tent, and the equation of s tate of the dark energy and, po- 
tentially, its evolution (e . g.. iHolder et al.ll200lt lHaiman et al.l 
[200lt IWeller et al.ll2002T: Maiumda r& Mohrll2003l) . In addi- 
tion, clusters probe the growth of structure in the Universe, 
which provides constraints different from and complemen- 
tary to the geometric con straints by the supernovae type la 
(e.g jAlbrecht et aLll2006l) . In particular, the constraints from 
structure growth may be crucial in distinguishing between the 
possibilities of the cosmic acceleration driven by dark energy 
or modification of the magnitude-redshift relation due to th e 
non-GR gravity on the largest scales (e.g.. lKnox et al.l l2005). 

The potential and importance of these constraints have mo- 
tivated current efforts to construct several large surveys of 
high-redshift clusters both using the ground-based optical 
and Sunyaev-Zel'dovich (SZ) surveys and X-ray missions in 
space. In order to realize the full statistical power of these sur- 
veys, we must be able to make accurate predictions for abun- 
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dance evolution as a function of cosmological parameters. 

Traditionally, analytic models for halo abundance as 
a function of m ass, have been used for estimating ex- 
pected evolution (iPress & Schechterlll974t iBond et al.lfl99ll 
iLee & Shandarinll998l:IShefh & Tormenl 19991) . Such models, 
while convenient to use, require calibration against cosmolog- 
ical simulations. In addition, they do not capture the entire 
complexity of halo formation and their ultimate accuracy is 
likely insufficient for precision cosmological constraints. A 
precision mass function can most directly be achieved through 
explicit cosmological simulation. 

The standard for precision deter mination of the mass func- 
tion from simulatio ns was set by iJenkins et al.l (1200 ll) and 
lEvrard et al.l d2002l) . who have presented fitting function for 
the halo abundance accurate to ~ 10-20%. These studies 
also showed that this function was universal, in the sense 
that the same function and parameters could be used to pre- 
dict halo abu n dance for different redshifts and cosmologies. 
IWarren et al.1 (H006) have further improved the calibration 
to w 5% accuracy for a fixed cosmology at z = 0. Several 
other studies have tested the universality of the mass function 
at high redshifts dReed et al.l l2003l 120071 iLukic et al. 1 120071: 
ICohn & Whitdl2007l) . 

One caveat to all these studies is that the theoretical counts 
as a function of mass have to be converted to the counts as 
a function of the cluster properties observable in a given sur- 
vey. Our understanding of physics that shapes these proper- 
ties is, however, not sufficiently complete to make reliable, 
robust predictions. The widely adopted strategy is therefore 
to calibrate abundance as a function of total halo mass and cal- 
ibrate the relation between mass and observable cluster prop- 
erties either sepa rately or within a survey itself using nuisance 
parameters (e.g. iMajumdar & Mohrll2004t iLima & Hull2004 
2005] |2003). The success of such a strategy, however, de- 
pends on how well cluster observables correlate with total 
cluster mass and whether evol ution of this corre lation with 
time is sufficiently simple (e.g jLima & Hull2005l) . 
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Tight intrinsic correlations between X-ray, SZ, and opti- 
cal observables and cl uster mass are expe cted theore tically 
(e.g.JBialek et al.ll200TtldaSilva et alJ200 4t lMotl et alJ200a 
iNagail 120061: iRravtsov et al.l 20061) and were shown to ex- 
ist observationally (e . g.. iMohr et alJ 119991 : iLin et al J |2004t 
I Vikhlinin et aH 120061: Maugh an 2007; Ar naud et alJ 120071: 
Isheldon et alJl2007l:IZhang et alj|2008h in the case when both 
observables and masses are defined within a certain spheri- 
cal radius enclosing a given overdensity. The majority of the 
mass function calibration studies, however, have calibrated 
the mass function with halos and masses measured using the 
friends-of-friends (FOF) percolation algorithm. This algo- 
rithm is computationally efficient, straightforward to imple- 
ment, and is therefore appealing computationally. The rela- 
tion between the FOF masses and observables is, however, 
quite uncertain. 

As we show below (see § 12.31 and Fig. [2j, there is large, 
redshift-dependent, and asymmetric scatter between the FOF 
mass and mass measured within a spherical overdensity, 
which implies that there is also large asymmetric scatter be- 
tween the FOF mass and cluster observables. This does not 
bode well for self-calibration of such relations. Furthermore, 
there is no way to measure the equivalent of the FOF mass 
in observations, which means that any calibration of the FOF 
mass and observables will have to rely on simulations. An 
additional issue is that halos identified with an FOF algorithm 
may not have one-to-one correspondence to the objects iden- 
tified in observational surveys. For example, the FOF finder 
is known to join neighboring halos into a single object even if 
their centers are located outside each others virial radii. Such 
objects, however, would be identified as separate systems in 
X-ray and SZ surveys. 

Although no halo-finding algorithm applied on simulations 
containing only dark matter may be perfect in identifying all 
systems that would be identified in a given observational sur- 
vey, the spherical overdensity (SO) halo finder, which identi- 
fies objects as spherical re gions enclosing a cer tain overden- 
sity around density peaks dLacev & Cold[l994h . has signifi- 
cant benefits over the FOF both theoretically and observation- 
ally. Most analytic halo models (see, e.g., ICooray & Shethl 
I2002L for review) assume that halos are spherical, and the 
statistics derived are sensitive to the exact halo definition. 
To be fully self-consistent, the formulae for halo properties, 
halo abundance, and halo bias, on which the calculations rely, 
should follow the same definition. The tight correlations be- 
tween observables and masses defined within spherical aper- 
tures means that connecting observed counts to theoretical 
halo abundances is relatively straightforward and robust. At 
the same time, the problem of matching halos to observed 
systems is considerably less acute for halos identified around 
density peaks, compared to halos identified with the FOF al- 
gorithm. 

Thus there is substantial need for a recalibration of the 
halo mass function based on the SO algorithm for a range 
of overdensities probed by observations and frequently used 
in theoretical calculations (~ 200 - 2000). Such calibra- 
tion for the standard ACDM cosmology is the main focus 
of this paper. Specifically, we focus on accurate calibra- 
tion of halo abundances for intermediate and high-mass halos 
(~ 10 11 - IO'^Mq) over the range of redshifts (z ~ 0-2) 
most relevant for the current and upcoming large cluster sur- 
veys. 

The paper is organized as follows. In § 2 we describe our 
simulation set and SO algorithm. In § 3 we present results for 



the mass function, demonstrating how our results depend on 
cosmology and redshift. In § 4 we summarize and discuss our 
results. 

Throughout this paper we use masses defined within radii 
enclosing a given overdensity with respect to the mean density 
of the Universe at the epoch of analysis. 

2. METHODS 
2.1. Simulation Set 

Table 1 lists our set of simulations. All the simulations are 
based on variants of the flat, ACDM cosmology. The cosmo- 
logical parameters for the majority of the simula tions reflect 
the Z eitgeist of the first-year WMAP results ( Sper geTet alJ 
120031) . We will refer to this cosmology as WMAP1 . A smaller 
number of simulations have cosmologi cal parameters clo ser 
to the three-year WMAP constraints JSpergel et alJl2007h . in 
which both Q m and cjg are lower and the initial power spec- 
trum contains significant tilt of n = 0.95. This subset of simu- 
lations are not of the same identical parameter set, but rather 
represent slight variations around a parameter set we will refer 
to globally as WMAP3. 

The largest simulation by volume followed a cubic box of 
1280 /i _1 Mpc size. There are fifty realizati ons of this simu - 
lation performed with the GADGET2 code (Springel 2005), 
which have been kindly provided to us by R. Scoccimarro. 
With the exception of these 1280 /i _1 Mpc boxes, the initial 
conditions for all simulations were created usi ng the stan- 
dard fi rst-order Zel'dovich approximation (ZA). ICrocce et al] 
(2006) point out possible systematic errors in the resulting 
mass function if first-order initial conditions are started in- 
sufficiently early. Using second order Lagrange perturbation 
theory (2LPT) to create initial conditions, they identify dis- 
crepancies betwee n the halo mass fun ction from their simu- 
la tions and that of IWarren et al.l d2006b at the highest masses. 
In lWarren et alJ d2006l) . several boxes larger than 768 h~ l Mpc 
were utilized in the analysis that are not listed in Table 1. 
In tests with our spherical overdensity halo finder, we find a 
discrepancy between the 2LPT simulations and these simula- 
tions. At this point, it is not yet clear whether the d iscrepancy 
is due to the effect advocated by lCrocce et al.l (l2006h or due to 
other numerical effects. We explore the issue of initial starting 
redshift in some detail in the Appendix A. What is clear, how- 
ever, is that results of these simulations systematically deviate 
from other higher resolution simulations, especially for larger 
values of overdensities. We therefore do not include them in 
our analyses. 

The first five simul ations listed in Table 1 were used 
in IWarren et al.l d2006l) in their analyses. The integrations 
were performed with the Hashed Oct-Tree (HOT) code of 
IWarren & Salmon! ([1993). Additionally, there are two HOT 
simulations in the WMAP3 parameter set. These simulations 
will be referred to in the text by their box size, in h~ l Mpc, 
prefixed by the letter 'H' . Simulations in the WMAP3 set will 
be appended with the letter 'W'. Due to identical box sizes 
between parameter sets, H384 will refer to the WMAP1 sim- 
ulation, H384W will refer to the simulation with WMAP3 
parameters, and H384f2 will refer to the low-il m simulation 
(which we will include in the WMAP3 simulation subset). 

There are six simu lations using th e Adap tive Refinement 
Technique (ART) of iKravtsov et all (119971) . and four that 
use GADGET2 in addition to the L1280 realizations. The 
L80 and L120 ART boxes modeling the W MAP1 cosmol- 
ogy are described in IKravtsov et al.l d2004f) and L250 simu- 
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FIG. 1. — A graphical key for the list of simulations in Table 1. The upper panel shows point styles for all the WMAP1 simulations ordered by the 
box size. Each simulation is represented with a different color, while different point types represent different numerical codes: circles=HOT, squares=ART, 
triangles=GADGET2. The lower panel plots all WMAP3 simulations, as well as H384f2, the low-f2,„ simulation. See Table 1 for the details of each simulation. 



TABLE 1 

Properties of the Simulation Set 



Lbox h 1 Mpc 


Name 


e h 1 kpc 




m p h 1 Mq 


(U m ,n b ,crg,h,n) 


Code 


Zi 


Zout 


^max 


768 


H768 


25 


1024 3 


3.51 X 10 10 


(0.3,0.04,0.9,0.7,1) 


HOT 


40 





800 


384 


H384 


14 


1024 3 


4.39 X 10' 


(0.3,0.04,0.9,0.7,1) 


HOT 


48 





3200 


271 


H271 


10 


1024 3 


1.54 X 10 9 


(0.3,0.04,0.9,0.7,1) 


HOT 


51 





3200 


192 


H192 


4.9 


1024 3 


5.89 x 10 8 


(0.3,0.04,0.9,0.7,1) 


HOT 


54 





3200 


96 


H96 


1.4 


1024 3 


6.86 x 10 7 


(0.3,0.04,0.9,0.7,1) 


HOT 


65 





3200 


1280 


L1280 


120 


640 3 


5.99 x 10" 


(0.27,0.04.0.9,0.7,1) 


GADGET2 


49 


0,0.5, 1.0 


600 


500 


L500 


15 


1024 3 x2 


8.24 x 10 9 


(0.3,0.045,0.9,0.7,1) 


GADGET2 


40 


0, 0.5, 1.25, 2.5 


3200 


250 


L250 


7.6 


512 3 


9.69 x 10' 


(0.3,0.04,0.9,0.7,1) 


ART 


49 


0,0.5, 1.25, 2.5 


3200 


120 


L120 


1.8 


512 3 


1.07 x 10 9 


(0.3,0.04,0.9,0.7,1) 


ART 


49 


0,0.5, 1.25, 2.5 


3200 


80 


L80 


1.2 


512 3 


3.18 x 10 s 


(0.3,0.04,0.9,0.7,1) 


ART 


49 


0, 0.5, 1.25, 2.5 


3200 


1000 


L1000W 


30 


1024 3 


6.98 x 10 1U 


(0.27,0.044,0.79,0.7,0.95) 


ART 


60 


0,0.5, 1.25, 2.5 


3200 


500 


L500Wa 


15 


512 3 x2 


6.20 x 10 10 


(0.24,0.042,0.75,0.73,0.95) 


GADGET2 


40 





3200 


500 


L500Wb 


15 


512 3 x2 


6.20 x 10 10 


(0.24, 0.042, 0.75, 0.73, 0.95) 


GADGET2 


40 





3200 


500 


L500Wc 


15 


512 3 x2 


6.20 x 10 10 


(0.24,0.042,0.8,0.73,0.95) 


GADGET2 


40 





3200 


384 


H384W 


14 


1024 3 


3.80 x 10 9 


(0.26,0.044,0.75,0.71,0.94) 


HOT 


35 





3200 


384 


H384f7„, 


14 


1024 3 


2.92 x 10 9 


(0.2,0.04,0.9,0.7,1) 


HOT 


42 





3200 


120 


L120W 


0.9 


1024 3 


1.21 x 10 8 


(0.27,0.044,0.79,0.7,0.95) 


ART 


100 


1.25,2.5 


3200 


80 


L80W 


1.2 


512 3 


2.44 x 10 8 


(0.23,0.04,0.75,0.73,0.95) 


ART 


49 


0, 0.5, 1.25, 2.5 


3200 



NOTE. — The top set of 5 simulations are from the Warren et al. 12006) study. The second list of 5 simulations are of the same WMAP1 cosmology, 
but with different numerical codes. The third list of 8 simulations are of alternate cosmologies, focusing on the WMAP3 parameter set. The HOT code 
employs Plummer softening, while GADGET employs spline softening. The force resolution of the ART code is based on the size of the grid cell at the 
highest level of refinement. A max is the highest overdensity for which the mass function can measured directly. Above this A, halo mass are inferred 
from the rescaling procedure in §2.3. A graphical key of this table is shown in FigurefT] 



lation is described by Tasitsiomi et al. (2008, in preparation), 
while the three WMAP3 boxe s are presented here . The L500 
simulations are de scribed in iGottlober & Yepesl (120071) and 
lYepes et all (12007I) 9 . These simulations contain equal num- 
bers of dark matter and SPH gas particles (without cooling). 
The ART and GADGET2 simulations will be referred by their 
box size with prefix 'L'. WMAP3 simulations have a 'W as 

9 see also http://astro.ft.uam.es/marenostrum/universe/index.html 



a suffix. 

Our simulation set comprises three different N-body codes, 
one based on the popular tree algorithm (HOT), one based 
on grid codes with small-scale refinement of high-density re- 
gions (ART), and one that combines grid and tree algorithms 
(GADGET2). We present a key in Figure [T] that graphically 
displays the range of box sizes. Each simulation is repre- 
sented by a different color, while different point types refer to 
different simulation codes: circles for HOT, squares for ART, 



4 



and triangles for GADGET2. These point symbols and colors 
will be used consistently in the figures below. 

2.2. Halo Identification 

The stand ard spherical overdens ity algorithm is described 
in detail in lLacev & Cole (1994). However, in our ap- 
proach we have mad e several important modifications. In 
lLacev & Cold (1994) the centers of halos are located on the 
center of mass of the particles within the sphere. Due to sub- 
structure, this center may be displaced from the main peak 
in the density field. Observational techniques such as X-ray 
cluster identification locate the center of the halo at the peak 
of the X-ray flux (and therefore the peak of density of the hot 
intraclustergas). Optical cluster searches will often locate the 
cluster center at the location of the brightest member, which 
is als o expected to be located near the peak of X-ray emis - 
sion (iLin et al.ll200H iKoester etaf]l2007l: iRykoff et al.ll2008h . 
Thus we locate the centers of halos at their density peaks. 

Our halo finder begins by estimating the local density 
around each particle within a fixed top-hat aperture with ra- 
dius approximately three times the force softening of each 
simulation. Beginning with the highest-density particle, a 
sphere is grown around the particle until the mean interior 
density is equal to the input value of A, where A is the 
overdensity within a sphere of radius R<\ with respect to 
the mean density of the Universe at the epoch of analysis, 
p m (z) = fi m (z)Pciit(z) = A»(0)(1 +z) 3 : 



(4/3)nR 3 A p m ' 

All values of A listed in this paper are with respect to p m (z). 

Since local densities smoothed with a top-hat kernel are 
somewhat noisy, we refine the location of the peak of the 
halo density with an iterative procedure. Starting with a ra- 
dius of r = /?a/3, the center of mass of the halo is calculated 
iteratively until convergence. The value of r is reduced iter- 
atively by 1% and the new center of mass found, until a fi- 
nal smoothing radius of R&/15, or until only 20 particles are 
found within the smoothing radius. At this small aperture, the 
center of mass corresponds well to the highest density peak of 
the halo. This process is computationally efficient and elim- 
inates noise and accounts for the possibility that the chosen 
initial halo location resides at the center of a large substruc- 
ture; in the latter case, the halo center will wander toward the 
larger mass and eventually settle on its center. Once the new 
halo center is determined, the sphere is regrown and the mass 
is determined. 

All particles within /?a are marked as members of a halo 
and skipped when encountered in the loop over all parti- 
cle densities. Particles located just outside of a halo can be 
chosen as candidate centers for other halos, but the iterative 
halo-centering procedure will wander into the parent halo. 
Whenever two halos have centers that are within the larger 
halo's Ra, the halo with the largest maximum circular veloc- 
ity, defined as the maximum of the circular velocity profile, 
V c (r) = [GM(< r)/r\ x l 2 , is taken to be the parent halo and the 
other halo is discarded. 

We allow halos to overlap. As long as the halo center does 
not reside within Ra of another halo, the algorithm identifies 
these objects as distinct structures. This is in accord with X- 
ray or SZ observations which would identify and count such 
objects as separate systems. The overlapping volume may 
contain particles. Rather than attempt to determine which 



halo each particle belongs to, or to divide each particle be- 
tween the halos, the mass is double-counted. No solution 
to this problem is ideal, but we find that the total amount of 
double-counted mass is only ~ 0.75% of all the mass located 
within halos, with no dependence on halo mass. This paral- 
lels the treatment of close pairs of clusters detected observa- 
tionally. When two X-ray clusters are found to have overlap- 
ping isophotal contours, each system is treated individually 
and double counting of mass will occur as well. 

For each value of A, the halo finder is run independently. 
Halo mass functions are binned in bins of width 0.1 in logM 
with no smoothing. Errors on each mass function are obtained 
by the jackknife method; each simulation is divided into oc- 
tants and the error on each mass bin is obtained through the 
variance of the halo number counts as each octant is rem oved 
from the full simulation volume (cf. IZehavi et al.1 12005, equa- 
tion [6]). The jackknife errors provide a robust estimate of 
both the cosmic variance, which dominates at low masses, 
and the Poisson nois e that dominates at high masses (see 
iHu & Kravtsovll2b03l for a the relative contributions of each 
source of error as a function of halo mass). 

When fitting the data, we only use data points with error 
bars less than 25% to reduce noise in the fitting process. We 
note that mass bins will be correlated (low-mass bins more 
so than high mass ones). We do not calculate the full covari- 
ance matrix of each mass function, so the x 2 values obtained 
from the fitting procedure should be taken as a general guide 
of goodness of fit, but not as an accurate statistical measure. 
However, we note that the data from multiple simulations in 
each mass range will be uncorrelated, and the lack of a covari- 
ance matrix should not bias our best-fit values for the mass 
function. 

2.3. Comparison ofFOF and SO halos 

IWhitd d2001l 12002b demonstrated that there is scatter be- 
tween the masses of halos identified with the FOF and SO 
halo definitions, as well as an offset between the mean halo 
masses using the canonical values of the linking length I = 0.2 
in the FOF algorithm and overdensity A = 200 in the SO ap- 
proach. Figure|2]compares the masses of halos identified with 
these two definitions for three different simulations. Halos in 
a simulation are first identified with our SO approach, then 
the FOF finder is subsequently run, beginning at the center 
of the SO halo. Figure^ compares A = 200 to / = 0.2. The 
symbols represent the median mass ratio tm = M200 /A^fof.2 as 
a function of M^oo. The curves represent the upper and lower 
90% bounds on the distribution of mass ratios. Although the 
median is near unity, the scatter is large and highly asymmet- 
ric. 

The asymmetry in the distribution is due to the FOF algo- 
rithm linking two or more distinct objects in close proximity 
to each other. Because we allow halos to overlap, FOF will 
treat these halos as a single object. Due to the arbitrary shape 
of FOF halos, the algorithm also links SO objects that do not 
overlap. The median mass ratio is also sensitive to the number 
of particles per halo; FOF halos are know n to be biased towar d 
higher masses at low particle number (Wa rren et alj 2006). 
The scatter between mass definitions is not alleviated by mak- 
ing the linking length smaller. This is shown in Figured, in 
which the same results are shown for A = 1600 and I = 0.1. 
The median is once again near unity, and the scatter remains 
identical. We note also that there is an offset in the median be- 
tween simulations as well; the results from L1000W are ~ 5% 
lower than the other simulations at I = 0.2 and ~ 10% lower 
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FIG. 2. — Comparison between spherical overdensity masses and friends-of-friends masses for the same sample of objects from H384, L250, and L1000W. 
Panel (a) compares the masses of A = 200 halos to FOF halos with / = 0.2. The symbols represent the median mass ratio, for objects binned by Miqo- The curves 
show the upper 90% and lower 10% bounds of the distribution of mass ratios in each A/200 bin: solid for H384, dashed for L250, and dotted for L1000W. The 
asymmetry in the mass ratio distribution reflects the tendency of FOF to link objects together. Panel (b) compares A = 1600 halos to FOF objects with 2 = 0.1. 
Panel (c) shows the distribution of mass ratios, vm =M2oo/^fof.2. for halos 13 < log M200 < 14 (solid line). The long tail of the distribution at ym < 0.5 indicates 
SO halos that are linked with other virialized objects in the FOF halo-finding process. The dotted line is the same distribution at z = 1.25. Panel (d) shows the 
distribution of i~m for the same mass range, for the A = 1600 and and FOF linking length / = 0.1. Solid and dotted lines are z = and z = 1.25, respectively. Both 
panels (c) and (d) show results for the L250 run. 



at / = 0.1. This offset is not due to the change in cosmology 
between the L1000W and the other simulations, therefore it 
must be a result of the lower mass resolution. 

We find that the curvature in the median mass ratio is alle- 
viated when adjusting the masses Mfof.2 by the Warren et. al. 
correction formula, (1 — A'" 6 ), where N p is the number of par- 
ticles in a halo. However, the curvature is not entirely amelio- 
rated by this formula at I = 0.1, demonstrating that the mass 
errors in FOF halos depend on the linking length. We find 
that (1— Np ) is sufficient to remove the FOF bias for / = 0.1. 
Figures |2l; and[2ji show the distribution of mass ratios for ha- 
los between 10 13 and 10 14 /! _1 M Q for one of the simulations. 
The solid histograms present results at z = and the dotted 
histograms is for z = 1 .25. Both the z = histograms exhibit a 



large, constant tail to low ratios. At higher redshift, the asym- 
metry of P{rM) becomes even stronger. The correlation be- 
tween spherically-defined masses and the FOF masses is thus 
broad and evolves with redshift. 

This has significant implications for comparisons with ob- 
servational cluster counts. Given that cluster observables 
correlate strongly with the spherical overdensity masses, the 
large scatter between Ma and Mfof indicates that the FOF 
correlation will be weaker. If one is to use a halo mass func- 
tion calibrated against halos and masses identified with the 
FOF algorithm, a significant additional effort would be re- 
quired to calibrate the scatter between FOF masses and ob- 
servables as a function of redshift, mass, and cosmology. In 
addition, this calibration will have to rely solely on theoretical 
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FIG. 3. — The halo density profiles are compared to analytic predictions 
for three different simulations. In each panel, the dotted curve repr esents the 
mean interior density given by an NFW profile with c(M) from D olag et all 
12004). The shaded region is the expected scatter assuming a\ ogc = 0.12. The 
solid curves with errorbars represent the numerical results. The left panel 
shows results from H384 for all halos M > 10 14 ' 5 h~ l Mq. The center and 
right panels show results for halos M > 10 15 h Mq. The center and left 
panel demonstrate that halo profiles are well resolved in these simulations. 
The right panel, shows significant deviations from the expected NFW profile 
at r < 0. l/?200 m the lower-resolution L1280 simulation. 
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Fig. 4. — Test of the resolution of the large-volume simulations, L500, 
L1000W, and one realization of LI 280. In each panel, the mass functions 
are plotted as residuals with respect to the best-fit f(cr) function from Ta- 
ble 2. The symbols represent the mass functions measured directly from 
the simulations at each A. The curves are mass functions inferred from the 
A = 200 halo catalog of each simulation, where the mass of each A = 200 
halo is scaled to higher overdensities assuming an analytic NFW halo (in- 
cluding scatter in concentrations at fixed mass). For the two higher resolution 
simulations, the scaled and true mass function are in agreement. Due to in- 
sufficient resolution, the L1280 mass function falls below the scaled mass 
function at high A. 



modeling, because the mass equivalent to the FOF cannot be 
directly measured in observations. The use of the halo abun- 
dance predictions made with the spherical overdensity algo- 
rithm is therefore strongly preferred. 

2.4. Accounting for effects of resolution 

Defining the halo masses within a radius enclosing a given 
overdensity stipulates that the halo mass is the integrated den- 
sity profile within a fixed radius. This means that the mass 
depends on the internal density distribution of the halo, and 
is thus more susceptible to the effects of resolution. The 
FOF masses, on the other hand, are measured within a given 
isodensity surface, and are therefore less sensitive to the in- 
ternal mass distribution. For example, iLukic et ail (2007) 
demonstrate that a reasonable FOF mass function can be ob- 
tained through a low-resolution simulation with as little as 8 
timesteps. If the same simulation is performed twice with dif- 
ferent resolutions, the same density peak in the lower resolu- 
tion simulation will have a shallower density profile and will 
in general have a different measured mass, Ma- The result is a 
systematic artificial shift in the estimated halo mass function. 
This effect will be larger for larger values of A, as smaller 
radii that enclose larger overdensities are more affected by 
resolution issues. 

To measure the SO mass function reliably at high A, we 
test whether the halo density profiles are properly resolved in 
each of the analyzed simulations at the overdensity in ques- 
tion. Figure [3] illustrates one of the resolution tests that we 
performed. It compares the halo density profiles from simu- 
lation s to the expected pro files. For the latter we use the well- 
tested [Navirroetal] d 1997b profile (hereafter NFW) with the 
concentration f or a given mass m easured in high-resolution 
simulations by Dola g et ail ((2004) 10 and a scatter in concen- 
tration of 0.12 in log 10 c. In this figure we show examples of 
one HOT simulation (H384), one ART simulation (L1000W), 
and one GADGET2 simulation (L1280). The HOT and ART 
simulations have force resolutions of 14 and 30 h~ l kpc, re- 
spectively, which is well within the scale radius of a typi- 
cal cluster-sized halo. The results for both the mean pro- 
file and its scatter are in excellent agreement with the NFW 
profile. The L1280 simulation has a force resolution of 120 
h~ l kpc, and deviations from the expected profile become clear 
at r < 0.1/?2oo. These differences will propagate into the esti- 
mate of the mass function if they are not taken into account. 

The results of comparisons similar to those shown in Figure 
[3] clearly identify which radii and which simulations profiles 
are affected by resolution. These results can then be used to 
determine the range of overdensities for which masses can be 
measured reliably in a given simulation. This is illustrated in 
FigureHJ which shows the mass functions from three different 
simulations at four values of A. The mass functions are plot- 
ted relative to the best-fit mass functions at each A, which are 
described in more detail below in § 3. At each overdensity we 
compare the mass functions measured in simulations to mass 
functions obtained by taking the individual halos found using 
the SO halo finder with A = 200 and rescaling their masses 
assuming the NFW p rofile, taking into account scatter in con- 
centrations (see, e.g.. lWhitdl200ll:lHu & Kravtsovll2003l) . We 
use the concentration-mass relation and scatter measured di- 

10 C200OT = 9.59 X (M/10 14 )-° 102 , normalized to the WMAP1 cosmol- 
ogy. When changing cosmology, we shift t he normalization using the frac- 
tional change in concentration from the Bullock et al. (2001) model at M = 
10 13 /T 1 M m . 
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rectly from our simulations (Tinker et. al., in preparation). 
The figure shows that the measured and re-scaled mass func- 
tions are in good agreement for A < 800, where the scaled- 
up mass function is ~ 5% higher than the true mass func- 
tion. This error is accrued from the halos located within R 200, 
which can become separate halos for higher overdensities and 
are not accounted for in the rescaling process. 

At higher overdensities, the agreement is markedly worse, 
especially for the lower-resolution L1280 boxes. At A = 
1600, the measured mass function is underestimated by ~ 
10%, increasing to ~ 20% at A = 3200. Therefore, for this 
simulation we use the directly-measured mass function only 
at A < 600, while at higher A we calculate the mass function 
by mass re-scaling using halos identified with an overdensity 
A = 600. A scaling baseline of log( Ahi g h / Ai ow ) < 0.9 accrues 
only < 2% error in the amplitude of the mass function at these 
masses. Thus the rescaled halo catalogs are reliable for cali- 
brating the halo mass function at high overdensity. This pro- 
cedure is used to measure high-A mass functions for L768 
(for A > 800) and L1280 (for A > 600). 

At A = 200 we choose a conservative minimum value of no 
less than 400 particles per halo. Below this value resolution 
effects become apparent, and simulations with differing mass 
resolutions begin to diverge . This is readily see n in the SO 
mass functions analyzed in iJenkins et al.l d2001l) . At higher 
A, halos are probed at significantly smaller radii, and the res- 
olution requirements are more stringent. Thus at higher A 
we increase the minimum number of particles such that, at 
A = 3200, A^ m i n is higher by a factor of 4. Exact values for 
each overdensity are listed in Table 2. 

3. HALO MASS FUNCTION 

3.1. Fitting Formula and General Results 

Although the number density of collapsed halos of a given 
mass depends sensitively on the shape and amplitude of the 
power spectrum, successful analytical ansatzes predict the 
halo abundance quite accurately by using a universal func- 
tion describing the mass fraction of matter in peaks of a given 
height, v = 5 c /a(M,z), in the linear density field smoothed 
at some scale R = (3M /4-KQ m ) 1/3 dPress & Schechtedll974t 
iBond et all [19911 ISheth & Tormenl 119991) . Here. S, ^ 1.69 
is a constant corresponding to the critical linear overdensity 
for collapse and a(M,z) is the rms variance of the linear den- 
sity field smoothed on scale R(M). The traditional nonlinear 
mass scale M* corresponds to a = S c . This fact has moti- 
vated the search for accurate universal f u nctions desc ribing 
simulation results b y IJenkins et al.l (1200 ll) . IWhitd (120021) . and 
IWarren et al.l d2006l) . Following these studies, we choose the 
following functional form to describe halo abundance in our 
simulations: 

dn p m dln<r~ l 



-0.64 



-0.52 



log(l/a) 
-0.38 -0.21 



-0.01 0.24 



0.55 



■ = /(f) - 

dM M dM 



(2) 



Here, the function f(a) is expected to be universal to the 
changes in redshift and cosmology and is parameterized as 



(- 

Kb 



+ 1 



-c/a 1 



where 



a = / P(k)W(kR)k z dk, 



(3) 
(4) 



and P(k) is the linear matter power spectrum as a function of 
wavenumber k, and W is the Fourier transform of the real- 
space top-hat window function of radius R. It is convenient to 
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FlG. 5. — The measured mass functions for all WMAP1 simulations, 
plotted as (M 1 /p m ) dn/dM against logM. The solid curves are the best-fit 
functions from Table 2. The three sets of points show results for A = 200, 
800, and 3200 (from top to bottom). To provide a rough scaling between 
M and <7~' , the top axis of the plot shows <r~' for this mass range for the 
WMAP1 cosmology. The slight offset between the L1280 results and the 
solid curves is due to the slightly lower value of Q„, = 0.27. 



recall that the matter variance monotonically decreases with 
increasing smoothing scale, thus higher M corresponds to 
lower a. 

The functional form (O was used in Warr en et al.l ((2006), 
with mi nor algebraic diff e rence , and i s similar to th e form s 
used by ISheth & Tormenl (119991) " and IJenkins etal.1 (120011) . 
Parameters A, a, b, and c are constants to be calibrated by 
simulations. The parameter A sets the overall amplitude of 
the mass function, while a and b set the slope and amplitude 
of the low-mass power law, respectively. The parameter c 
determines the cutoff scale at which the abundance of halos 
exponentially decreases. 

The best fit values of these parameters were determined 
by fitting eq. (f3]l to all the z = simulations using x 2 mini- 
mization and are listed in Table 2 for each value of A. For 
A > 1600, we fix the value of A to be 0.26 without any loss of 
accuracy 12 . This allows the other parameters to vary mono- 
tonically with A, allowing for smooth interpolation between 
values of A. 

Figure|5]shows the mass function measured for three values 
of A and the corresponding best fit analytic functions. We plot 
(M 2 /p,„) dn/dM rather than dn/dM to reduce the dynamic 
range of the y-axis, as dn/dM values span nearly 14 orders 
of magnitude. The figure shows that as A increases the halo 

1 1 A convenient property of the Sheth & Tormen mass function is that 
one recovers the mean matter density of the universe when integrating over 
all a~ l . Equation (5J does not converge when integrating to <r~' = 0. In 
Appendix C we present a modified fitting function that is properly normalized 
at all A but still produces accurate results at z = 0. 

12 Although a four-parameter function is required to accurately fit the data 
at low A, at high overdensities the error bars are sufficiently large that a 
degeneracy between A and a emerges, and the data can be fit with only three 
free parameters, given a reasonable choice for A. 
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FIG. 6. — Panel (a): The measured f(cr) from all simulations in Table 1. Results are presented at z = and for A = 200. The solid line is the best fit function of 
equation (5). The lower window shows the percentage residuals with respect to the fitting function. In the WMAP1 cosmology, the range on the data points on 
the jc-axis is roughly 10 105 h~ l Mq to 10 15 ' 5 hT 1 Mq. Panel (b): The measured f(a) at z = 1.25. We restrict results to simulations for which we have previous 
redshift outputs. In the WMAP1 cosmology, the range of data points on the ;t-axis is 10" hT 1 M© to 10 15 hT 1 Mq. The solid line is the same as in panel (a), 
which was calibrated at z = 0. The lower window shows that the z = 1.25 mass function is offset by ~ 20% with respect to the results at z = 0. 



masses become systematically smaller. Thus from A = 200 
to 3200, the mass scale of the exponential cutoff reduces sub- 
stantially. The shape of the mass function is also altered; at 
A = 200 the logarithmic slope at low masses is ~ -1 .85, while 
at A = 3200 the slope is nearly -2. This change in slope is due 
to two effects. First, the change in halo mass accrued with 
changing the halo definition A is not independent of mass. 
Because halo concentrations depend on mass, c/Mai does not 
equal dM^i for halos of two different masses. 

Second, a number of low-mass objects within /?200 of a 
larger halo are "exposed" as distinct halos when halos are 
identified with A = 3200. Although all halos contain sub- 
structure, these "revealed" subhalos will only impact overall 
abundance of objects at low mass, M< 10 12 /t'Mq, because 
the satellite fraction (the fraction of all halos located within 
virial radii of larger halos) decre ases rapidly from » 20% to 
zero for M > 10 12 /t'Mq (e.g. iKravtsov et al.ll2004 . This 
trend can be understood using average properties of subha- 
los in parent CDM halos. Subhalo populations are approx- 
imately self-similar with only a we ak trend with mass (e.g., 
iMoore et all 1 19991: iGao et al.fl2004l) . and the largest subhalo 
typically has a mass of rs 5- 10% of the host mass. Thus, 



at a given mass M only hosts with masses > 10M can pro- 
duce significant number of new halos when halo identifica- 
tion at higher A is performed. At high masses, the number of 
such halos decreases exponentially with mass, and therefore 
the contribution of such "exposed" halos becomes small. 

Figure [6^ shows the function f(a) measured for all simu- 
lations in Table 1 at z = with A = 200. The solid curve is 
equation (fJJ using the best-fit parameters from Table 2. The 
residuals with respect to this fit demonstrate the high accu- 
racy of our numerical results and the consistency of different 
codes, mass resolutions, and cosmologies. Figure [6J) shows 
f(a) at z = 1 .25 for a subset of simulations for which higher 
redshift outputs are available. The solid curve represents the 
results from z = 0. At this redshift, the results at ~ 20% below 
the z = results, nearly independent of a~ l . This demonstrates 
that the mass function is not universal in redshift, or for cor- 
respondingly large changes in cosmology, 13 at this level of 
accuracy. We address evolution of f(a) with z in §3.3 below. 

13 Note that we can interpret higher redshift outputs of a given simulation 
as the z = epoch of a simulation with different cosmological parameters 
corresponding to f2 m (z) and other parameters at the redshift in question. 
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FIG. 7. — Residuals of the measured mass functions with respect to the best fit analytic mass functions from Table 2 for all WMAP1 simulations at z = 0. Error 
bars are shown for the first and last point for each simulation, and o nly points with less than 10% error bars are plotted, with the exception of L80, for which 15% 
is the maximum. For A = 200, the blue curve represents the Jenkins et al. 12001) SO180 mass function (scaling up to A = 200 yields indistinguishable results). 
The red dashed curve represents the Sheth & Tormen 1 1999) mass function. For A = 400, the blue curve represents the Jenkins et al. (2001) S0324 (scaled up to 
A = 400). For A = 1600 and A = 800, the solid curve represents the Jenkins SO(324) mass function scaled up analytically assuming NFW profiles. 



TABLE 2 

Mass Function Parameters for f(a) at z = 



A 


A 


a 


b 


c 


X 2 /u (ALL) 




X 2 /v (WMAP1) 


X 2 /v (WMAP3) 


X 2 /v (WMAP3-fit) 


200 


0.186 


1.47 


2.57 


1.19 


1.15 


400 


1.07 


1.66 


1.62 


300 


0.200 


1.52 


2.25 


1.27 


1.17 


400 


1.08 


1.65 


1.60 


400 


0.212 


1.56 


2.05 


1.34 


1.05 


600 


0.96 


1.49 


1.37 


600 


0.218 


1.61 


1.87 


1.45 


1.06 


600 


0.99 


1.55 


1.28 


800 


0.248 


1.87 


1.59 


1.58 


1.10 


1000 


1.07 


1.36 


1.14 


1200 


0.255 


2.13 


1.51 


1.80 


1.00 


1000 


0.97 


1.22 


1.16 


1600 


0.260 


2.30 


1.46 


1.97 


1.07 


1600 


1.03 


1.34 


1.25 


2400 


0.260 


2.53 


1.44 


2.24 


1.11 


1600 


1.07 


1.50 


1.26 


3200 


0.260 


2.66 


1.41 


2.44 


1.14 


1600 


1.09 


1.61 


1.33 



NOTE. — A^m is the minimum number of particles per halo used in the fit. Fits are for simulations at z = 0. 
The WMAP1 and WMAP3 X 2 l v values are with respect to the WMAP1 and WMAP3 simulations, respectively, but 
using the best-fit parameters. The WMAP3-fit X 2 l v values are independent fits using only the WMAP3 simulations 



3.2. Results as a function of A 

The best-fit parameters of equation (f3]l resulting from fits to 
all simulations for 9 values of overdensity are listed in Table 



2. Figure [7] shows the residuals of individual WMAP1 simu- 
lations with respect to global fits at different A. We include 
L1000W in these plots to show consistency between cosmolo- 
gies at cluster masses. For the fifty realizations of L1280, we 
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FIG. 8. — The trajectories of the best-fit parameters of f(a) from Table 2. In each panel, the order of the points is from low-A to high-A (left to right). 
Error bars represent l-cr variance of parameters from the MCMC chain. In panel (a), the normalization A is plotted against power-law slope a. In panel (b), the 
power-law amplitude b is plotted against the cutoff scale c. The lower panel shows the rms scatter of mass functions from 100 bootstrap samples, creating by 
sampling the simulation list. Light gray is for A = 200, while dark gray is for A = 1600. 



plot the mean /(er) and the error in the mean. Each panel 
shows the fractional residuals of the measured mass functions 
with respect to the best-fit /(er) for four values of A. To avoid 
crowding, error bars are plotted for the maximum and mini- 
mum mass scale for every simulation; the latter is represen- 
tative of the cosmic variance given the finite simulation vol- 
ume, while the former is dominated by Poisson noise. We 
list formal values of \ 2 l v for ° ur diagonal error bars in Ta- 
ble 2. The values in column 6 are for all z = simulations, 
while the value in column 8 is the \ 2 l v f° r the same param- 
eters but with respect to the WMAP1 simulations only. Not 
surprisingly, the \ 2 l v values reduce slightly when compar- 
ing the best-fit f(a) to the WMAP1 simulations only, which 
comprise a majority of the simulations and therefore drive the 
fitting results. 

The solid blue curve in the A = 200 panel represents the fit- 
ting function of [Jenkins et al.l (120011) calibrated on their set 
of rCDM simulations (their equation B3), using A = 180 
(rescaling this equation to 200 yields nearly indistinguishable 
results). At M > 10 12 fi "'M fl , the Jenkins resul t is 10-15% 
below our results. The Sheth & Tormen (1999) function is 
similarly offset from our results. In the A = 400 panel, the 
blue curve shows the Jenkins et. al. fitting function calibrated 
to A = 324 on their set of ACDM simulations (essentially 
the WMAP1 cosmology). For this comparison the Jenkins 



formula has been rescaled to A = 400 usin g the same halo 
rescali ng techniques discussed in §2.3 and in iHu & Kravtsovl 
(2003). The Jenkins SO(324) function (their Equation B4) 
is in good agreement with our results for M < 10 13 h~ l M , 
while at higher masses there are variations of ±5-10%. 

The solid curves in the A = 800 and 1600 panels are 
the Jenkins SO(324) result scaled up to those overdensities. 
At cr" 1 > 1, the residuals increase, while at lower masses 
the rescaled f(a) underestimates the numerical results by 
5 - 10%. Both of these effects are due to subhalos becoming 
exposed when halos are identified using higher overdensity. If 
a high-mass halo contains a large subhalo, the rescaling pro- 
cedure will overestimate the mass of that object at higher A. 
At low masses, the rescaling procedure does not account for 
the revealed substructures. The change in mass from A = 200 
to A = 1600 is - 50% at 10 14 /i- 1 M . If subhalos are dis- 
tributed within parent halos in a similar fashion to the dark 
matter, then the rescaling procedure should underestimate the 
mass function by ~ 0.5 x 0.2 = 0. 1 (where 0.2 is the s ubhalo 
fraction for low-mass halos from lKravtsov et al. 2004). 

Figure[8]shows that the best fit parameters of f(a) vary with 
A smoothly. This means that interpolating between these 
best-fit parameters can be expected to yield accurate mass 
function parameters at any desired overdensity. In Appendix 
B we show examples of the interpolated mass functions, as 
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FIG. 9. — Residuals of the measured mass functions with respect to the best fit analytic mass function from Table 2 for all WMAP3 simulations at z = 0. Error 
bars are shown for the first and last point for each simulation, and only points with less than 10% error bars are plotted, with the exception of L80W, for which 
15% is the maximum. 



well as fitting function for the f(a) parameters as a function 
of A. The error bars are la and are obtained by marginalizing 
over all other parameters. The errors on the amplitude A are 
i~ 3 - 4%, but this parameter is highly correlated with b and 
the true scatter about the best-fit /(er) is < 1 % at most masses. 

The lower panel in Figure [8] shows the rms scatter in our 
constraints on f(a). The scatter was calculated by bootstrap 
resampling of the simulation set and repeating the fitting pro- 
cess on 100 bootstrap samples. 14 The shaded area is the vari- 
ance of the bootstrap fits. The light gray region represents 
results for A = 200, while the dark gray region represents 
A = 1600. Between er -1 = 0.63 and er -1 = 1 .6 the scatter is less 
than 1% (M= 10 11 5 h~ l M Q and 10 15 r'M s for the WMAP1 
cosmology). Outside this mass range the results diverge due 
to lack of coverage by the simulations. Because the WMAP1 
simulations dominate by number, these constraints should be 

14 Because the 50 realizations of L1280 outnumber all the rest of the sim- 
ulations (which only number 17), we create bootstrap samples by first sam- 
pling from the list of LI 280 realizations, then sampling from the rest of the 
simulation set. This guarantees a fair sampling of the range of <t~' probed 
by the simulations. If we do not do this, many bootstrap samples will only 
contain mass function with results above M > 2 X 10 14 h Mq, which would 
artificially inflate the size of the low-mass errors. 



formally regarded as the accuracy of the fit for the WMAP1 
cosmology. 

Figure|9]compares the calibrated mass functions from Table 
2 with the measured mass functions from the WMAP3 sim- 
ulations (i.e., the last seven entries in Table 1). Column 9 of 
Table 2 contains values of x 2 l v f° r the WMAP3 simulations 
only. The X 2 l v ^ somewhat larger than for the WMAP1 
runs at all overdensities, even though the WMAP3 residuals 
do not seem to be systematically offset from the global /(er) 
fits. We test this statistically by refitting for the parameters of 
equation (0) using only the WMAP3 simulations. The x 2 l v 
values are listed in column 10 of Table 2. For each A, the 
X 2 of the fit is only reduced marginally. In the mass range 
covered by our simulations, the difference between the global 
f(a) functions and those derived from the WMAP3 simula- 
tions differ by < 2%, but with a ~ 4% uncertainty in the nor- 
malization of the WMAP3-only fitting function, derived from 
the bootstrap method described above. Thus we conclude that 
the higher x 2 values are not due to a systematic change in 
/(cr) due to variations in cosmology, but rather scatter in the 
simulations themselves at the ^5% level, excluding obvious 
outliers where Poisson noise dominates. 
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FIG. 10. — Redshift evolution of the A = 200 mass function. Each panel 
shows the residuals of the z = mass function with respect to the measured 
mass functions at z = 0, 0.5, 1.25, and 2.5. Note that the simulation set used 
here is a combination of WMAP1 and WMAP3 boxes. Error bars are shown 
for the first and last points for each simulation, and only points with < 10% 
are shown, with the exception of the L80 and L80W, for each 15% is the 
limit. The shaded region brackets 10 13 hT 1 M Q to 10 14 hT 1 Mq. The solid 
curves represent the z. = mass function modified by equations (5) — {8}. 
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FIG. 1 1 . — Redshift evolution of the A = 1600 mass function. Each panel 
shows the residuals of the z = mass function with respect to the measured 
mass functions at z = 0, 0.5, 1.25, and 2.5. Note that the simulation set used 
here is a combination of WMAP1 and WMAP3 boxes. Results are plotted 
down to halos with 400 particles, as opposed to the limit of 1600 used in 
fitting /(cr). All points with errors < 15% are plotted. The shaded region 
brackets 10 3 h" 1 Mq to 10 14 h" 1 Mq. The solid curves represent the z = 
mass function modified by equations (3} — (8). 



3.3. Redshift Evolution 

FigureQJJshows the evolution of the A = 200 mass function 
for four different redshifts from z = to 2.5. Results are plot- 
ted for the subset of simulations for which we have previous 
redshift outputs. When modeled as pure amplitude evolution, 
the mass function evolves as (1 +z)~ a26 . However, it is also 
clear that the shape of /(cr) is evolving with redshift such that 
the amplitude at a~ l > 1 decreases at a higher rate. This is 
more evident in Figure[TT] in which /(<j) at A = 1600 is shown 
for the same redshifts. As A increases, both the evolution in 
the amplitude and shape of /(cr) become stronger. 

In Figures [10] and QT| the solid curves show a model in 
which the first three parameters of /(cr) are allowed to vary 
as a power law of 1 +z; 



loga(A) 



A{z)=A a (\+z) 
a(z) = ao (l+zj 
b(z) = bo(l+zT l 
0.75 



-0.14 



0.06 



log. A 75)) 



(5) 
(6) 
(7) 

(8) 



where subscript '0' indicates the value obtained at z = in Ta- 
ble 2. Modulation of A controls the overall amplitude of /(er), 
while a controls the tilt, and b sets the mass scale at which the 
power law in f(a) becomes significant. Modifying b results 
in a shift between the amplitudes at low and high <r~ l , thus it 
encapsulates the changes in f(a) with A seen in Figures [10] 
and[TT] Although the redshift scaling introduced here matches 
the results at z < 2.5 accurately, residuals of > 5% emerge at 
z = 2.5. It is possible that the evolution between z = 1 .25 and 
2.5 is slowing down. Because the numerical results at z = 2.5 
are quite noisy and cover only a small range in cr" 1 , our results 
at this value of z and extrapolation to higher redshifts must be 
checked with other simulations. Extrapolating equation (0)- 
© to z = 10 produces an /(cr) that is reduced by ~ 50% with 
respect to z = 0. This seems unlikely given current studies but 
n eeds to be checked with a consistent halo finding algorithm. 

Reed et all d2007t) parameterize the redshift-dependent 
mass function in terms of both a and the effective spectral in- 
dex of the linear power spectrum, n e g. These authors use this 
parameterization to model the mass function at z > 10, where 
differences in the slope of n e ff from z = are large. This ap- 
proach is ill suited for modeling the evolution at z < 3, where 
there is very little change in the effective spectral index. 

It is interesting to note that the evolution in the exponen- 
tial cutoff scale is minimal. Any evolution in this mass scale 
would yield quantitatively different residuals than those seen 
in Figure [10] and [TT] Namely, the residuals would show pro- 
nounced curvature at cr" 1 > 1. Our results show that the dom- 
inant effect is a shift in the normalization in the mass function 
rather than the cutoff mass scale. Thus our results are not 
consistent with f(a) being universal as a function of virial 
overdensity because A v j, evolves with redshift. Nor are our 
results consistent with the mass function being universal at a 
fixed overdensity with respect to the critical density (rather 
than defining A with respect to the background, as we do 
here). 

The Jenkins et al. (2001) study reports no detected evolu- 
tion of the FOF or SO mass functions with redshift. More re- 
cent results quantify the evolution of the FOF at high redshift , 
z > 10, to be 5-10 % dLukic et all 120071: iReed et alJ 120071: 
Cohn & White 2007). However, friends-of-friends identified 
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FIG. 12. — Evolution of the FOF mass function for linking lengths of / = 0.2 (left panels) and 1 = 0.1 (right panels). The simulations used are L500 and L250. 
The L500 simulation has been downsampled to 1/8 the original particle number. For / = 0.2, the residuals are plotted with respect to the Warren et al] 120061) 
function. Mass functions are plotted down to 100 particles per halo but have not been corrected for discreteness effects (i.e., equation 2 in Warren et. al.). For 
/ = 0.1, the residuals are plotted with respect to the A = 1600 mass function from Table 1. Note the larger range of the y-axis at z = 2.5 for 2 = 0.1. The FOF mass 
function evolves less than the SO mass function, but this largely a numerical effect due to increased linking of distinct halos. 



halos may have a different response to changes in the redshift 
evolution of halo profiles. Merging rates vary with redshift, 
and this may be reflected in the FOF tendency to bridge dis- 
tinct structures. Figure [12] shows the redshift evolution for 
friends-of-friends selected halos in the L500 and L250 boxes. 
The panels in the left column show the results for halos iden- 
tified with a linking length of 0.2. Res iduals are calculated 
with respect to the Warr en et all (|2006) fitting formula with 
their best fit parameters, plotted down to halos containing 100 
particles. The friends-of-friends masses have not been cor- 
rected for any systematic errors (equation [2] in I Warren et all 
2006), resulting in the slight negative slope to the residuals at 
low masses. The mass function shows some redshift evolu- 
tion, but only of order ~ 10% at z= 1.25, or roughly half that 
in Figure [10] 

The right column shows the results for halos identified with 
a linking length of / = 0. 1 . The smaller linking length iden- 
tifies halos with higher overdensities. The residuals are with 
respect to f(a) for A = 1600. For this linking length, the red- 
shift evolution is stronger than for I = 0.2, and the shape of the 
FOF mass function changes dramatically. As a whole, these 
results indicate that the mass function is also non-universal 



for FOF halos, with the degree of non-universality depending 
on the linking length used. 

These results are in general agreement with those of other 
recent studies that considered evolution of the mass function 
for FOF halos, although the overall picture of how the mass 
function evolves with redshi ft is not yet clear. The simulation 
results of iLukic et all d2007l) exhibited ~ — 5% residuals with 
respect to the z = Warren et. al. mass function as z = 5, but 
with a monotonic trend of rising residuals with increasing red- 
shift. The FOF mass function in the Millennium Simulation, 
show s roughly 20% evolution from z = to 10 dReed et all 
12001) . Finally, iFakhouri & Mai (120071) recently showed that 
the Millennium Simulation FOF mass function, once cor- 
rected for spurious FOF linking between halos, evolves by 
~ 20% from z = to 1 . This is consistent with our findings, 
but note that the volume of the Millennium simulation and 
statistics at large masses is substantially worse than in our set 
of simulations. 

4. SUMMARY AND DISCUSSION 

We have presented a new fitting function for halo abun- 
dances and their evolution in the ACDM cosmology. The fit- 
ting function can be used to predict halo mass functions for 
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FIG. 13. — Halo mass range corresponding to the range of cr~' on which 
/((j) is calibrated. The shaded region bounded by the solid curves shows how 
this mass range evolves with redshift for the WMAP1 cosmology. The dashed 
curves show the upper and lower mass limits for the WMAP3 cosmology of 
the L80W simulation. The dotted line indicates the maximum redshift output 
of our simulation set. 



spherical aperture masses defined with an arbitrary overden- 
sity over a wide range of values. For the WMAP1 cosmology 
our results are accurate at the percent level in the mass range 
relevant for cluster cosmology. For the WMAP3 cosmology 
our results are accurate to < 5%. One of our main results is 
that the mass function is non-universal, and varies in a sys- 
tematic way with redshift in the interval z = [0,2.5], with the 
abundance of halos at a given cr" 1 monotonically decreasing 
with increasing z. 

We have parameterized redshift evolution of /(cr) as a sim- 
ple scaling of the z = fitting parameters with (1 +z) a - We 
note that if this evolution is driven by changes in f2 m with 
z, it may be more robust to model f(a,z) as a function of 
the growth rate rather than 1 + z. Our simulation set does not 
probe a large enough cosmological parameter space to detect 
differences due to different growth factors. However, this will 
become important when investigating how the mass function 
evolves in dark energy cosmologies, in which the primary 
change in structure formation is a different growth function 
of perturbations. 

Our finding of evolving, non-universal /(cr) is quantita- 
tively different from the results of previous analyses that use 
the friend-of-friends method for halo identification, which 
generally show weaker evolution and greater degree of univer- 
sality of the function /(cr). We argue that the likely explana- 
tion for this difference is greater sensitivity of the SO defined 
mass to the redshift evolution of halo concentrations. As dis- 
cussed previously, SO masses are the integrated halo profiles 
within a specified radius and lower halo concentrations result 
in lower masses at fixed abundance (or, conversely, fewer ha- 
los at fixed mass). The fact that the high-mass end of the 
mass function (where concentrations at z = are lower and 
the mass within /?20o/c20o is a significant fraction of the total 



mass) evolves somewhat faster than the low-mass end, argues 
that evolution of concentrations plays a significant role in the 
evolution of /(cr). 

The evolution of halo concentrations is mostly driven by 
the change in fi m with redshift. This implies that /(cr) in cos- 
mologies with substantially different matter densities at z = 
will be systematically different from the one we find here (per- 
haps closer to our z > 1 results). There are indications that this 
is indeed the case. The H384S1 sim ulation, with f l „ = 0. 2, is 
above /(cr) by - 5% at z = 0. The Uenkins et al] d2001h fit- 
ting function for A = 180 was calibrated on simulations with 
n m = 1, producing a fit ~ 15% below our results at the same 
overdensity. The Jenkins SO(180) mass function is close to 
our A = 200 results at z = 1 .25, where Q m is approaching unity. 

The lower evolution of the FOF mass function with redshift 
can be understood from Figure [2] The distribution of mass 
ratios between FOF and SO halos changes between z = and 
z= 1.25. The median mass ratio, Mso/A^fof, decreases while 
the scatter increases at higher z due to more linking of dis- 
tinct objects. The number of distinct objects at a fixed cr" 1 
decreases, but the higher incidence of linking offsets this ef- 
fect. Thus the weaker evolution of the FOF mass function is 
due to this linking of separate collapsed halos and is largely 
artificial. The better universality of /(cr) may still seem like 
an advantage of the FOF mass function. However, as we dis- 
cussed in this paper, the large and redshift-dependent scat- 
ter between SO and FOF masses implies similarly large and 
redshift-dependent scatter between FOF masses and cluster 
observables. This makes robust interpretation of observed 
cluster counts in terms of the FOF halo counts problematic. 

Our fitting function is calibrated over the range 0.25 < 
cr" 1 < 2.5, which at z = spans a range of halo masses roughly 
10 10 - 5 < M < 10 15 5 h~ x Mr, depending on the specific choice 
of cosmology. In Figure Qj] we show how this mass range 
evolves with redshift. By z = 3, the lower mass limit is ~ 10 5 
/i _1 M . At this redshift, ou r fitting function is in agreement 
with the numerical results of Colin et al.] (120041) . which probe 
the mass range 10 5 < M < 10 9 /T 1 M Q . At higher redshifts, 
cr" 1 is a slowly varying function of mass, making the lower 
mass limit evolve rapidly. Because our calibration of the 
redshift dependence of the mass function parameters extends 
only to z = 2.5, we caution against extrapolation of equations 
© — (© to significantly higher redshifts. As noted above, 
/(cr) is evolving less rapidly from 1.25 < z < 2.5 than from 
< z < 1 .25. Thus using the z = 2.5 /(cr) should yield a mass 
function with reasonable accuracy at higher z, but must be 
verified with additional simulations. 

The range of cosmologies probed here is narrow given the 
volume of parameter space, but it is wider than the allowed 
range given recent results from CMB in combi nation with 
other large-scale measures (Komatsu et al. 2008). For gen- 
eral use that does not require 5% accuracy, extending our re- 
sults somewhat outside this range will produce reasonable re- 
sults. It is unlikely that variations in the shape and amplitude 
of the power spectrum will yield significantly different forms 
of /(cr). As discussed above, however, large variations in O m 
at z = (ie, Q m = 0.1 or il m = 1), are not likely to be fit by our 
z = mass function within our 5% accuracy. Models with a 
higher matter density at z = can be approximated by using 
our calibrated /(cr) at the redshift for which f2 m (z) is equal to 
the chosen value. 

The next step in the theoretical calibration of the mass func- 
tion for precision cosmology should include careful exam- 
ination of subtle dependencies of mass function on cosmo- 
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logical parameters (especially on the dark energy equation of 
state), effect s of neutrinos with non-zero mass, eff ects of non- 
gaussianitv (iGrossi et al.ll2007t iDalal et al.ll2007l) . etc. Last, 
but not least, we need to understand the effects of baryonic 
physics on the mass distribution of halos and relate d effects on 
the m ass function, whi ch can be quite signi ficant ( Rud d et all 
120081) . The results of lZentner et alj d2007l) indicate that the 
main baryonic effects can be encapsulated in a simple change 
of halo concentrations, which would result in a uniform shift 
of Ma and a uniform correction to f(a). Whether this is cor- 
rect at the accuracy level required remains to be demonstrated 
with numerical simulations. 

Our study illustrates just how daunting is the task of cal- 
ibrating the mass function to the accuracy of < 5%. Large 
numbers of large-volume simulations are required to esti- 
mate the abundance of cluster-sized objects, but high dynamic 
range is required to properly resolve their internal mass dis- 
tribution and subhalos. The numerical and resolution effects 
should be carefully controlled, which requires stringent con- 
vergence tests. In addition, the abundance of halos on the ex- 
ponential cutoff of the mass function can be influenced by the 
choice of method to generate initial conditions and the start- 
ing re dshift, as was recently demonstrated by ICrocce et al.1 
(2006, see also Appendix A). All this makes exhaustive stud- 
ies of different effects and cosmological parameters using 
brute force calibration of the kind presented in this paper for 
the ACDM cosmology extremely demanding. Clever new 
ways need to be developed b oth in the choice o f the param- 
eter space to be investigated (Habib et al. 2007) and in com- 



plementary studies of various effects using smaller, targeted 
simulations. 
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FIG. 14. — Comparison between the large-box simulations used in the text and those in Warren et al. 12006) and Evrard et al. 12002). The box sizes and 
point-types for the three HOT boxes and the Hubble Volume are shown in the top panel above the horizontal line. In addition, a version of the L1000W ART box, 
started at lower redshift, is also included in the comparison. The large-box simulations used from Table 1 are also included below the horizontal line. The bottom 
panel compares the A = 200 mass functions to the best-fit /(it) at z = 0. The middle panel shows the results at z = 1.25. In the z = panel, the shaded region 
indicates I0 l4 h~ l M© < M < 10 15 h~ l Mq in the WMAP1 cosmology. In the z = 1.25 panel, the shaded region indicates 10 14 /r' Mq < M < 10 14 ' 5 Mq. 
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APPENDIX 

A. TESTS OF THE INITIAL CONDITIONS 

In a recent study. ICrocce et al] J2006) investigated differences between using the standard first-order Zel'dovich Approximation 
(ZA) and second-order Lagrangian perturbation theory (2LPT) for generating initial conditions of cosmological simulations. ZA 
assumes that particle trajectories are straight lines, but for large density fluctuations trajectories should curve due to tidal effects. 
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Thus, if a simulation is initialized at the epoch where the overdensity is lar ge in some region s, the resulting error in particle 
trajectories will lead to 'transients' in the evolution of perturbations (see also Scoccimar rol 19981) . which can persist to z = 0. This 
effect is strongest for the regions containing rarest peaks of largest height that tend to evolve into the largest galaxy clusters at low 
redshift. In their simulation results, Crocce et. al. find a ~ 10% discrepancy at M ~ 10 15 h~ l M in z = mass functions between 
2LPT and ZA with starting redshift of z; = 24. This discrepancy is expected to grow more significant at higher redshift at fixed 
halo mass. The effect is particularly worrisome for precision calibration of abundance of the most massive objects at any redshift 
(those objects that are currently collapsing or have only recently collapsed). In this appendix we present tests of the effects of the 
initial redshift on the mass function and explain why we have discarded some of the large-volume simulations from our analysis. 

The top panel in Figure [14] shows a graphical key of the three large-box HOT simulations used in the Warren et. al. fit that 
we do not utilize in our mass function fits. These simulations have starting redshifts of Zi = 34, 28, and 24 (with z; decreasing 
with increasing box size). In addition, we also have results from the Hubble Volume (HV ) simulation, a 3000 h Mpc simulation 
(lEvrard et al.l l 2002l We use the same SO halo catalog presented in lEvrard et"aTI (120021) . which used a density criterion of 200 
times the critical density rather than the mean. Thus we have scaled the halo masses from A = 666 to A = 200, assuming NFW 
profiles as detailed in §2. Lastly, we have included a re-simulation of the L1000W ART box which has been initialized at z, = 35 
rather than z; = 60 using the same set of random phases and ZA at both starting redshifts. 

The bottom panels of Figure [14] show the residuals of the simulation mass function from the best fit to our core simulation set 
at z = and z = 1.25. At 10 14 h M©, all simulations are in excellent agreement. However, at 10 15 h~ l M Q , the HOT boxes are 
~ 10-20% below the f(a) obtained from the 2LPT simulations and ART L1000W run. The mass function of the HV simulation, 
with Zi = 35, is also ~ 15% below the 2LPT simulations. 

At z = 0, there is a ~ 2% difference betwe en low-z; ART box and the higher-z, version used in the fitting. This is smaller than 
the difference between mass functions for the lCrocce et al.| [2006 simulations with z,- = 24 and z ( = 49, which may be due to sample 
variance. However, the difference between the two ART boxes increases at higher redshift. The ART box with z, = 60 is in good 
agreement with the 2LPT simulations at z = 1 .25, implying that convergence has been reached at a lower z ( than shown in Crocce 
et al. The run with z, = 35, however, is 20-40% lower than the best fit at large masses. 

It is not yet entirely clear whether the source of the discrepancies in the mass functions at the highest masses can be attributed 
solely to the errors of the ZA-generated initial conditions. The difference between the large-volume HOT boxes and the 2LPT 
results are larger than expected from just the ZA errors. Also, both ART boxes, with z, = 35 and z, = 60, are in agreement with 
the 2LPT simulations at z = 0. Other factors, such as resolution effects on the halo density profiles, may play a dominant role in 
the discrepancy exhibited by both the HOT boxes and the HV simulation. Regardless of the source of the discrepancy, it is clear 
that the large-volume HOT boxes and HV simulations are systematically different from other higher-resolution simulations. We 
therefore do not include them in our analyses. 

In summary, the simulations which we use to derive our constraints on the high-mass e nd of the halo mass fu nction are all robust 
against changing initial redshift. The 2LPT simulations have been thoroughly tested in lCrocce et al.1 (120061) . The L1000W and 
L500 simulations, utilizing ZA with z; > 50, show consistent results with the 2LPT simulations at multiple redshifts. However, 
quantifying the effects of initial conditions, finite simulation volume, and possible numerical artifacts at the < 1% level will 
require significant additional work. 

B. INTERPOLATION OF MASS FUNCTION PARAMETERS 

To facilitate the use of our results in analytic calculations, we provide fitting functions for the parameters of f(a) as a function 
of log A. The dependence of each parameter in the mass function is reasonably well described by 

. _ / 0.1 (log A) -0.05 if A < 1600 n 
A ~|0.26 if A > 1600, ( } 

a= 1.43 + (logA-2.3) L5 , (B2) 

b= l.O + aogA-1.6)" 1 ' 5 , (B3) 

and 

c= 1.2 + (logA-2.35) L6 . (B4) 

All logarithms are base 10. Because the parameters of f(a) are not completely smooth with log A, these functions yield mass 
functions that are accurate to only < 5% for most values of A, but can degrade to < 10% at a~ x > 0.2 for some overdensities. 
Figure IBTBI demonstrates the accuracy of the fitting functions with respect to the results from Table 2. For higher accuracy, we 
recommend spline interpolation of the parameters as a function of log A. Figure iBTBl shows the results of the spline interpolation 
when obtaining the para meters of f(tr). We provide the second derivatives of the f(a) parameters for calculation of the spline 
coefficients (cf., §3.3 in lPress et all 1992b in Table B3. 

C. AN ALTERNATE, NORMALIZED FITTING FUNCTION 

The fitting function given in equation (01 is an excellent descriptor of the data over the range of our data, but at a~ l < 0.1, f(a) 
asymptotes to a constant value. For some applications, specifically halo model calculations of dark matter clustering statistics, it 
is necessary to integrate over all a~ l to account for all the mass in the universe. Integrating (0 over all cr" 1 yields infinite mass. 
In this appendix we present an alternative fitting function that is normalized such that 
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FIG. B 15. — Accuracy of the fitting functions presented in Appendix B for calculating the parameters of /(cr) as a function of A (solid lines). All curves are 
residuals with respect to the best-fit results of /(cr) from Table 2. For all overdensities except A = 600, the accuracy of /(cr) is < 5%. The dashed lines show the 
accuracy of /(cr) when using spline interpolation, which is accurate to < 2% for all A and cr~' . 



TABLE B 3 
Second derivatives of f(a) parameters 



A A a be 



200 0.00 0.00 0.00 0.00 

300 0.50 1.19 -1.08 0.94 

400 -1.56 -6.34 12.61 -0.43 

600 3.05 21.36 -20.96 4.61 

800 -2.95 -10.95 24.08 0.01 

1200 1.07 2.59 -6.64 1.21 

1600 -0.71 -0.85 3.84 1.43 

2400 0.21 -2.07 -2.09 0.33 

3200 0.00 0.00 0.00 0.00 



Jg(a)d\na~ l = 1 (CI) 

for all values of A at z = 0. We focus on equation © for our main results because the parameters of that function vary more 
smoothly and monotonically with A, and incorporating redshift evolution into that function is more straightforward and more 
accurate. Because we can only calibrate our mass function to cr" 1 > 0.25, the behavior of the fitting function at lower masses is 
arbitrary. Thus it is not to be expected that the fitting function in this appendix is more or less accurate than equation (0 below 
this calibration limit, merely that the function is better behaved. 
With these caveats in mind, we find that at z = a function of the form 

(C2) 



8(<r) = B 



- +CT 
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TABLE C4 

Normalized Mass Function Parameters for 
g{a) at z = 



A 


B 


d 


e 


/ 


8 


X> 


200 


0.482 


1.97 


1.00 


0.51 


1.228 


1.14 


300 


0.466 


2.06 


0.99 


0.48 


1.310 


1.16 


400 


0.494 


2.30 


0.93 


0.48 


1.403 


1.04 


600 


0.494 


2.56 


0.93 


0.45 


1.553 


1.07 


800 


0.496 


2.83 


0.96 


0.44 


1.702 


1.09 


1200 


0.450 


2.92 


1.04 


0.40 


1.907 


1.00 


1600 


0.466 


3.29 


1.07 


0.40 


2.138 


1.07 


2400 


0.429 


3.37 


1.12 


0.36 


2.394 


1.12 


3200 


0.388 


3.30 


1.16 


0.33 


2.572 


1.14 



yields nearly identical results to tho se pr esented in Figure [7] Equation iC2\ has four free parameters, with B set by the nor- 
malization constraint from equation tCli , We follow the same procedure for fitting the model to the data as in §2.4. Best-fit 
parameters are listed in Table C4. The x 2 / l/ values are similar to the values listed in Table 2. The asymptotic slope of f(a) in 
the lSheth & Tormenl (1 19991) fitting function is cr~° 4 at low masses. The values of / in Table C4 vary around this value, with the 
A = 200 g(a) going as cr~ ' 51 and A = 3200 going as cr" 033 . 

Another requirement of the halo model is that dark matter be unbiased with respect to itself. This requires a recalibration of 
the large-scale halo bias function, which we investigate in another paper (Tinker et al., in preparation). 



